A GENERALIZED FINITE ELEMENT METHOD 
FOR THE DISPLACEMENT OBSTACLE PROBLEM OF 
CLAMPED KIRCHHOFF PLATES 



SUSANNE C. BRENNER, CHRISTOPHER B. DAVIS, AND LI-YENG SUNG 

Abstract. A generalized finite element method for the displacement obstacle problem of 
clamped Kirchhoff plates is considered in this paper. We derive optimal error estimates and 
present numerical results that illustrate the performance of the method. 



1. Introduction 

Let fl be a bounded polygonal domain Q C M 2 , / G L 2 (fl), g G H 4 (Q), and ipi,ip2 £ 
C 2 (Q) PI C(Cl) be two obstacle functions such that 

(1.1) ipi < ifj 2 in f2 and ipi < g < 1P2 on dfl. 
Consider the following problem: Find u G H 2 (Q) such that 

(1.2) u = argminG(v), 

where 

(1.3) K = {ve H 2 {Vl) : u-ge H 2 (n) } ^i < v < ip 2 on SI}, 

(1.4) G(v) = ^a(v,v)-(f,v), 

(1.5) /'VV:VW, f/.,) /'/,,/, 

and V 2 f : V 2 w = Ylij=i v XiXjW XiXj is the (Frobenius) inner product of the Hessian matrices 
of v and w. 

Since K is a nonempty closed convex subset of H 2 (Q) and a(-, •) is symmetric and coercive 
on Hq(Q) which contains the set K — K, it follows from the standard theory [281 EH EH [22] 
that (1.2) has a unique solution u G K characterized by the following variational inequality: 

(1.6) a(u, v — u) > (f, v — u) Wv<eK. 
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The convergence of finite element methods for second order obstacle problems was inves- 
tigated in [THi [131 E] j shortly after it was shown in [TT] that the solutions for such obstacle 
problems belong to H 2 (Q) under appropriate regularity assumptions on the data. This full 
elliptic regularity allows the complementarity form of the variational inequality (in the strong 
sense) to be used in the convergence analysis. 

In contrast, it was shown in EEE EH EH] that the solution u of ( |1 .2 ) / ( 1.6 ) belongs to 



Hf oc (Q) fl C 2 (Q) under the assumptions above on /, g, if)j and ip2- Since the obstacles are 



separated from each other and from the displacement boundary condition (cf. ( 1.1 )), we have 
A 2 u = f near dQ. Therefore it follows from the elliptic regularity theory for the biharmonic 
operator on polygonal domains [HI EH [El EI] that u G H 2+a (N) for some a G (|, 1] in an 
open neighborhood M of dVt. The elliptic regularity index a is determined by the interior 
angles of Vt and we can take a to be 1 for convex Vt. Thus the solution u of ( 1.2[ )/( jL6 ) 
belongs to H 2+a (Q) fl Hf^Q) fl C 2 (Q) in general. Moreover, it is easy to construct examples 
where u Hf oc (p) even for smooth data [TH] 



This lack of Hf oc (p) regularity means that the complementarity form of (1.6) only exists 
in a weak sense [15J. Consequently convergence analysis based on the weak complementarity 



form of (1.6) would only lead to suboptimal error estimates. 



A new convergence analysis for finite element methods for ( 1.2 )/( 1.6 ) that does not rely 
on the complementarity form of the variational inequality (1.6) was proposed in pi)], where 



optimal convergence was established for C 1 finite element methods, classical nonconforming 
finite element methods, and C° interior penalty methods for clamped plates (g = 0) on convex 
domains. The results in [10] were subsequently extended to general polygonal domains and 
general Dirichlet boundary conditions for a quadratic C° interior penalty method [9] and a 
Morley finite element method [5] . The goal of this paper is to extend the results in [HI E] to 
a generalized finite element method for plates [17, 31J. 

The rest of the paper is organized as follows. We introduce the generalized finite element 
method in Section [2] and carry out the convergence analysis in Section [3j Numerical results 
are reported in Section |4| 

2. A Generalized Finite Element Method 



We begin with the construction of the approximation space Vh in Section 2.1 and define 
an interpolation operator from H 2 (Q) into Vh in Section 2.2 The discrete obstacle problem 



is given in Section 2.3 We refer the readers to [31 E] for various aspects of generalized finite 
element methods. 

2.1. Construction of the approximation space. The approximation space is based on 
partition of unity by flat-top functions j29j [32] . 

2.1.1. Partition of Unity. Let be the C 1 piecewise polynomial function given by 



[x 



lR, 



[X) 

x) 



(l + x) 2 (l 
(1 -x) 2 (l 



2x) 
2x) 



if x G [-1,0] 
if x G [0, 1] 
if \x\ > 1 
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which enjoys the partition of unity property that 

(2.1) L (x - 1) + (f) R {x) = 1 for < x < 1. 

We define a flat-top function ip$ by 

' <P L ( X ' ( 2S +5) ) Hxe [-1-6,-1 + 5} 



ips{x) 



I^R ( X-{l-S) 

25 



ifze [-1 + 5,1-5] 

if x e [1 - 5, 1 + 5] 

if x £ [-1-5,1 + 5] 

Here 5 is a small number that controls the width of the flat-top part of this function where 
^ = 1. 

For ease of presentation we take Q to be a rectangle (a, b) x (c,d). But the construction 
and analysis can be extended to other domains (cf. Remark 2.3 and Examples 4 and 5 in 
Section fib. 



We first expand to a larger rectangle f2 = (a — 71, 5 + 71) x (c — 72, d + 72) where 
71 and 72 are two positive numbers, and then we divide Q into disjoint congruent closed 
rectangular patches Qj (cf. Figure 2.1) with center yj = (3/7,1,2/7,2), width hi and height h 2 , 
for j = 1, . . . , N . We assume that the numbers 

8j = 7j/{hj/2) (j = 1,2) 

belong to the interval [/3i,/3 2 ], where /3i and /3 2 are constants that satisfy < (3i < (3 2 < 1- 
For each patch Qj, let 



x\ - Vj,i 
hi/2 



1p6 



X2 ~ Vj,2 

h 2 /2 



It follows from (2.1) that 



1, . . . , N} is a partition of unity in Q, i.e., 

N 

y^^j = l on 0,. 

3=1 



The flat-top region of each patch, defined by 

Qf* = {a; G : ^(x) = 1}, 

is the rectangle centered at yj with width hi(l — 5i) = hi—2j 1 and height ^(l - 6 2 ) = h 2 —2^ 2 
(cf. Figure 2.1 ). 



Remark 2.1. By construction we have (cf. Figure 2.1) 

• Qj at n Qf at = if i ^ j. 

• The support of extends a horizontal distance of 71 = 5\{h/2) and a vertical 
distance of 72 = 5 2 {h/2) outside of the patch Qj. Hence the supports for and tyj 
will intersect in a rectangular region of width 271 or 272 if Qi is a neighbor of Qj. 

• If Qj n dVl ^ 0, then Qf 1 n dfi ^ 0. 
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Figure 2.1. Partitioning of a square domain f2: First we expand the shaded 
region f2 to f2 (left). Next we divide (l into congruent rectangles Qj, 1 < j ' < 16 
(middle). The flat top regions Q^ at , 1 < j < 16 are shown as the darker shaded 
regions (right). 



2.1.2. Approximation space. The space Q 2 of biquadratic polynomials will serve as the local 
approximation space and the global approximation space is defined to be 

N 

j'=i 

Below we present an explicit basis of Vh that will be used in our numerical computations. 
On the reference interval [—1,1] we have two types of quadratic polynomials: 

• Lagrange interpolation polynomials that satisfy N^(Lj) = Sij for 1 < i,j ' < 3, 
where JVi(u) = v(-l), N 2 (v) = v(0), and N 3 (v) = v(l). 

• Hermite interpolation polynomials Hi(£) that satisfy Ni(Hj) = 5ij for 1 < i,j < 3, 
where Ni(y) = v'(-l), N 2 (v) = v(-l), and N 3 (v) = v{l). 

The tensor product of different combinations of these polynomials will provide local bases 
on the two-dimensional rectangular patches. 
Let Tj : R 2 — > R 2 be defined by 

= (yj,i + Uh/2)(1 - 5 1 ),y 3 , 2 + 6(V2)(i - S 2 )). 

Then Tj maps the reference square [—1,1] x [—1,1] to the flat-top region Q^ at . 

Depending on the location of the patch Qj, we use different reference basis functions. 
There are three possibilities. 

• For those patches such that Qj D d£l = 0, the reference basis functions are 

(2.2) = i = 3(Z-l) + fc, l<A;,/<3. 

• For those patches such that intersects the boundary on only one side, say the 
vertical edge x± — a of f2, the reference basis functions are 

(2.3) ^(0 = Hki&Lfo), i — 3(1-1) + k, 1 < k, I < 3. 
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Note that in this case Tj maps the line £1 = —1 to the part of Qj that intersects d£l. 
The cases where Qj intersects other sides of fl can be treated analogously. 

• For those patches such that Qj intersects a corner of Q, say the lower left corner 
(a, c), the reference basis functions are 

(2.4) f Jt (0 = H^Hfo), i = 3(1 — 1) + k, l<k,l<3. 

Note that in this case Tj maps the corner (—1,-1) of the reference square to the 
lower left corner (a, c) of Q. The cases where Qj intersects other corners of Q can be 
treated analogously. 

The nodal variables (or degrees of freedom) for the local approximation space are depicted 



in Figure 2.2, where pointwise evaluations of functions, directional derivatives, gradients 
and mixed second order derivatives are represented by solid dots, arrows, circles and double 
arrows respectively. 




Figure 2.2. Reference elements described in (2.2), (2.3), and (2.4) respectively. 



An explicit basis for the global approximation space Vh is then given by 
(2.5) [j^oTr 1 ] : j = l,2,...,iV;z = l,2,...,9}. 

Remark 2.2. Since all the nodes in a rectangular patch Qj are located in Qj at where ^ = 
for k 7^ j, all the basis functions of V/, vanish at the nodes in Qj except those associated 
with Qj. 



Figure |2.3| illustrates the degrees of freedom associated with the basis of the global ap- 
proximation space for a square which is divided into 9 square patches where hi = h 2 = h 
and 5i = 5 2 = S. 

Remark 2.3. One may follow the same procedure for non-convex polygonal domains. As an 



example, consider an L-shaped domain 



a, a) 



2 \[0,a] 2 . One could divide the domain into 



rectangular patches everywhere except near the reentrant corner. Near the reentrant corner, 
one could construct local biquadratic polynomial basis functions in the reference L-shaped 
domain (— 1, 1) 2 \[0, l] 2 dual to the nodal variables 
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Figure 2.3. A partition of a domain Q differing only by the choice of S. The 
solid lines separate the different patches Qj,j = 1,...,9. The dashed lines 
represent the extension of Qj by 5(h/2) on each side. This figure also shows 
the location of the degrees of freedom corresponding to 5 = 1/6 (left) and 
5 = 1/3 (right). 



N 1 (v)=v(0,0), N 2 (v) 
N A (v)=v(-l,-l), N 5 (v) 

N 7 (v) = ^(0,0), N 8 (v) 



depicted in Figure |2.4 



v(l,0), N 3 (v) 
dv 

Wi 

dv 
W2 



-(0,0), N 6 {v) 
-(0,1), N 9 (v) 



v(0,l), 
dv 

Wi 

d 2 v 



(0,1), 
(0,0), 



Figure 2.4. Reference element for the L-shaped domain. 
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2.2. Interpolation Operator. First we define interpolation operators associated with the 
rectangular patches. Let £ G H 2 (M. 2 ). 



For a patch with the local basis given by (2.2) (cf. the reference element on the left 
of Figure 2.2), we define Uj( to be the polynomial in Q 2 such that (Hj£) °Tj = C °Tj 
at the 9 points in the set {(p, q) : p, q — —1, 0, 1}. 



For a patch with the local basis given by (2.3 ) (cf. the reference element in the middle 
of Figure 2.2), we define Uj( to be the polynomial in Q 2 such that 

(i) (rLyC) oTj = (oTj at the 6 points in the set {(p, q) : p = —1,1 and q = —1, 0, 1}. 

(ii) The polynomial (d[(Ilj() o Tj]/<9£i)|£ 1= _i equals the quadratic polynomial A(£ 2 ), 
which is the L 2 projection of (d(( o 2y)/^£0l6=-i into the space of quadratic 
polynomials in the variable £ 2 . 



• For a patch with the local basis given by (2.4) (cf. the reference element on the right 
of Figure 2.2), we define Uj( to be the polynomial in Q 2 such that 

(i) (HjC) ° Tj = ( o Tj at the 4 points in the set {(p, q) : p, q = ±1}. 

(ii) (d[(Uj() o Tj)/d£i)\£ 1= -i equals the quadratic polynomial A at £ 2 = ±1, where 
A(£ 2 ) is the L 2 projection of (d(( o Tj)/9^i)|^ 1= _i into the space of quadratic 
polynomials in the variable £ 2 . 

(iii) (d[(Jlj() o Tj]/d£, 2 )\&=-i equals the quadratic polynomial p, at £i = ±1, where 
/t(£i) is the L 2 projection of o I})/9^ 2 )|^ 2 =-i the space of quadratic 
polynomials in the variable £i. 

(iv) The value of (<9 2 [(IL,C) o T^/d^d^) at (-1,-1) equals (A'(-l) + p!{-\))/2. 

Remark 2.4. Since Tj maps the reference square to Qj, the interpolant Uj( is determined 
by the restriction of ( to Q^ at . 

We can now define the global interpolation operator LT^ : H 2 (Q) — > Vh by 

N 

n h ( = X;(n y Ct)*i vc€i/ 2 (fi), 

i=i 

where Cf G P 2 (M 2 ) is any extension of £. The interpolant 11^ is independent of the choice of 
£t by Remark |2.4[ Moreover, by construction we have 



(2.6) u h (H 2 (n)) = v h nH 2 {n). 

Let Qj be the rectangle centered at %jj with width h\{l + <5i) = /ii + 271 and height 
/i 2 (l + 5 2 ) = h>2 + 27 2 . Let /i = max(/ii, /i 2 ). Since ELjP = P for any P e Q 2 , the estimate 

2 

m=0 



cS 
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follows from the Bramble-Hilbert lemma [7J [18] and scaling. From here on we use C to 
denote a generic positive constant that is independent of the mesh size h. 



Combining the local interpolation error estimate (2.7) and the estimates for the partition 



of unity functions tyj in [32], we immediately have (cf. [291 132] ) the following error estimates 
for the global interpolation operator IT^: 



(2.8) 



m=0 



Remark 2.5. Classical rectangular C l finite element methods would require a local approxi- 
mation space that is at least bi-cubic [6] . Of course we can also use bi-cubic polynomials as 
the local approximation space in our GFEM (cf. jTTl ETj and Example 1 in Section [1]). 

2.3. The Discrete Obstacle Problem. Let Vh be the set of the nodes in the rectangular 
patches corresponding to the degrees of freedom involving pointwise evaluation of the local 



basis functions. (Such nodes are represented by solid dots in Figure 2.3 and Figure 2.4 ) 
The GFEM for the model problem is to find Uh G Kh such that 



(2.9) 



Uh 



argminG(f ) 

v£K h 



where the quadratic functional G is defined by (1.4)— (1.5) and 

(2.10) K h = {veV h :v-n h geH*(n),Mp)<v(p)<Mp) VpGV h }. 



Remark 2.6. Approximation of the essential boundary conditions u = g and du/dn = dg/dn 
are both included in the definition of Kh- Moreover Kh is nonempty because HhK C Kh by 



(1.3) and (2.6). 



Remark 2.7. In view of Remark |2.2| and the defining properties of the polynomials and 
Hi, the constraints defining Kh are box constraints with respect to the basis of Vh defined 



in (2.5). 



It follows from the standard theory that the discrete obstacle problem (2.9) has a unique 
solution characterized by the discrete variational inequality 

(2.11) a(u h ,v - u h ) > (f,v - Uh) Vv e K h - 



3. Convergence Analysis 



We begin with some preliminary estimates in Section 3.1 and introduce an auxiliary ob- 



stacle problem in Section 3.2 that connects the continuous problem (1.2) and the discrete 



problem (2.9). The main result is derived in Section 3.3 



3.1. Preliminary Estimates. In view of (2.8), it suffices to find an optimal estimate for 
\HhU — Uh\H 2 {n)- Using the discrete variational inequality (2.11), we have 

|ILyu - u h \ 2 H 2 {n) = a(U h u - u, ILvu - u h ) + a(u - u h , ILyu - u h ) 

< \Il h u - u\ H 2 {n) \Il h u - u h \ H 2(n) + a(u,U h u - u h ) - (f,Yl h u-u h ) 
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1 



U\ 



\ii 2 (n) + ^l^hU - u h \ 2 H 2(n) + a(u,U h u - u h ) - (f,U h u - u h ), 
which implies 

(3.1) \U h u - u h \ 2 H2{n) < \TL h u - u\ 2 H2{n) + 2 [a(u, U h u - u h ) - (f, U h u - u h )] . 

We can therefore complete the error analysis by finding an optimal estimate for the expression 
a(u, U h u - u h ) - (/, U h u - u h ). 



The following result is useful for the error analysis in Section 3.3 



Lemma 3.1. There exists a positive constant C independent of h such that 



(3.2) \a{<f>X-Il h ()\ <Ch 2a \\<f>\\ H 2 +a{Q) 

for all <j) G H 2+a (tt) and ( G H 2+a (n) n H 2 (tt). 

Proof. Let ( G H 2+a (VL) nH 2 (tt) be arbitrary. On the one hand we have an obvious estimate 



(3.3) 



|a(0,C-IIfcC)| < |0|fla(n)|C-n fc C|fla(n) 
< Ch a \4>\ H 2(n)\(\H 2 +<*{n) 



V0 G H 2 (Q) 

that follows from (|2.8|). On the other hand, we have another estimate 

(V ■ V 2 ^ 



|a(0,c-n fc c)| 



v(c - n h Q dx 



(3.4) 



< C\4>\ H 3 {n) \( - IlftClfli(n) 
<L7/i 1+Q |0|^ (n) |C|^(n) 



V0 G H 3 (Q) 



that follows from (2.6), (2.8) and integration by parts. 

The estimate (3.2) follows from (3.3), (3.4) and interpolation between Sobolev spaces 

HIES!- □ 

3.2. An Auxiliary Obstacle Problem. We can connect the continuous obstacle problem 
( 1.2 ) and the discrete obstacle problem (2.9 ) through an intermediate obstacle problem: Find 
Uh G Kh such that 

(3.5) Uh = argminG(v) 

veK h 

where 

(3.6) K h = {ve H 2 {9) :v-ge H 2 ^),^) < v(p) < Mp) G V h }. 

Note that Kh is a closed convex subset of H 2 (Q) and K C Kh- The unique solution of 
( 3.5[ ) is characterized by the variational inequality: 

(3.7) a(u h ,v-u h )>{f,v-u h ) V«e4 



The connection between (1.2) and (3.5) is given by the following properties of Uh from 



(3i 



\u - Uh\H 2 (n) 



< Ch, 
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and there exists h > such that 

(3.9) u h = u h + 5 h ,\(j)i - S h>2 (p2 £ K Vh<h , 

where <pi and 02 are C°° functions with compact supports in fl such that 0« = 1 on the 
coincidence set {iGll: u(x) = ^i(x)}, and the positive numbers 5h,i and 5h,2 satisfy 

(3.10) S h)i < Ch 2 . 



Note that v-U h ge H 2 (Q) for all v G K h (cf. ( |27lOp ) and hence, by Q , 

(3.11) u + - n h (y) G ^ A Vvelif fc . 

3.3. Error Estimates for the Generalized Finite Element Method. We now complete 
the error analysis of the generalized finite element method by deriving an optimal estimate 
for the expression a(u, Uh — Uh) — (/, n^n — Uh). To simplify the presentation, we introduce 
the transitive relation A < B defined by 

A<B^A-B< C(h 2a + h a \Tl h u - u h \ H 2 (u) ). 

Since 

a(u, Tl h u - u h ) = a(u - u h , H^w - u h ) + a(u h , TlhU - u h ) 

and 

a(u - u h ,U h u - u h ) <\u- W/i|i?2(n)|n ft M - u h \ H ^(u) < Ch\U h u - Uh\m{Q) 
by the estimate (3.8), we have 

(3.12) a(u, U h u - u h ) - (/, U h u - u h ) < a(u h , U h u - u h ) - (/, U h u - u h ). 



In view of (3.11), we can use the auxiliary variational inequality (3.7) to obtain 

Uh 



a(u h ,U. h u 



= a(u h , u h -u h -(g- Il h g)) + a(u h} U h u - u h + (g - U h g)) 
< (/, u h -u h - (g- Tl h g)) + a(u h , U h u -u h + {g- Uhg)), 



which together with (3.12) implies 
(3.13) a(u,U h u-u h ) - (f,U h u-u h ) ^ 

a(u h , Yl h u -Uh + (g- Khg)) - (/, w - Uh) - (/, (IlftW -u) + (g- U h g)). 



We can rewrite the first term on the right-hand side of (3.13) as 
a(u h , U h u -tih + {g- n h gf)) = a(u h - u, U h u -u h + {g- Uhg)) + a(u, U h (u - g) - (u- g)) 

+ a(u, u — Uh). 

Observe that 

a(u h - u, U h u -u h + (g- U h g)) = a(u h - u, (U h u - u) + (u - u h ) + (g - Tlhg)) 

< \uh ~ u\ H 2 {n) (\U h u - u\ H 2( n ) + \u- Uh\H 2 (Q) + \g~ n^|// 2 (n)) 

< Ch l+a 
by dZ8b and (Q, and 



a(u,U h (u-g) 



u 



g)) < Ch 



2a 
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by Lemma 3.1. Moreover we have, by (2.8), 

|£2(n)(||n fc «-«||i a (n) 



9 ~ n^|| La (n)) < Ch 
a(u,u- u h ) - (f,u- u h ). 



2+a 



-(/, (U h u-u) + (g-U h g)) < 

Combining these relations and (3.13), we arrive at the estimate 

(3.14) a(u, U h u - u h ) - {f, n h M - u h ). $ 

According to ( 1.6[ ), (3.9) and (3.10), we have 

a(u,u- u h ) - (f,u- u h ) 

= a(u, u-u h ) - (f, u-u h ) + S h ,i[a(u, 0i) - (/, 4>i)} - S h)2 [a(u, 2 ) - (/, ^2)] 
< Ch\ 



and hence (3.14) leads to the estimate 

a(u, U h u - u h ) - (/, H-hU - u h ) $ 0, 

which means 

(3.15) a(u, U h u - u h ) - (f, U h u - u h ) < C(h 2a + h a \U h u - u h \ HH 

Theorem 3.2. There exists a positive constant C independent of h such that 

\u - u h \ H 2 (n) < Ch a . 



i 2 (ft), 



Proof. It follows from (2.8), (3.1), (3.15) and the arithmetic and geometric means inequality 
that 



\K h u - u h \ 2 mm < C(h 2a + h a \U h u - u h \ H 2 (n) ) < Ch 2a + l\U h u - u h \ 2 m 



□ 



which implies 

(3.16) \U h u-u h \ H 2 (n) <Ch a . 
The theorem follows from ( 2.8[ ), (3.16) and the triangle inequality. 

Since H 2 (Q) is embedded in C(Q) by the Sobolev embedding theorem [HE3], the following 
corollary is immediate. But numerical results in Section [4] indicate that the convergence rate 
in the L 00(H) norm should be higher than the convergence rate in the H 2 (Q) norm. 

Corollary 3.3. There exists a positive constant C independent of h such that 



(3.17) 



\u - u h \ Loo(n) < Ch c 



Remark 3.4. Under additional assumptions 



on the exact coincidence sets (resp. free 



boundaries), the error estimate (3.17) implies the convergence of the discrete coincidence 
sets (resp. free boundaries) to the exact coincidence sets (resp. free boundaries). Details 
can be found in [9]. 



12 



SUSANNE C. BRENNER, CHRISTOPHER B. DAVIS, AND LI-YENG SUNG 



4. Numerical Results 

We present numerical results for several one-obstacle problems to demonstrate the perfor- 
mance of the GFEM. The obstacle function from below will be denoted by ip. The first four 
examples are from [9] . The discrete obstacle problems are solved by a primal dual active set 
strategy from [U [25] . 

Example 1. Here we apply the GFEM to a problem with a known exact solution to 
validate the numerical results. We begin with the plate obstacle problem on the disc {x : 
|x| < 2} with / = 0,ip(x) = 1 — \x\ 2 and homogeneous Dirichlet boundary conditions. This 
problem is rotationally invariant and can be solved exactly. The exact solution is 



u(x) 



Ci\x\ 2 In \x\ + C 2 \x\ 2 + C 3 ln \x\ + C 4 , r < \x\ < 2 



1 — |x| 2 , \x\ < r 



where r w 0.18134452, C x « 0.52504063, C 2 « -0.62860904, C 3 « 0.01726640, and C 4 
1.04674630. We then consider the obstacle problem on Q = (—0.5, 0.5) 2 whose exact solution 
is the restriction of u to Q. For this problem / = 0, ip(x) = 1 — \x\ 2 and the (non-homogeneous) 
Dirichlet boundary data are determined by u. 



We partition Cl following the procedure described in Section \2.1\ and define j to be the 
level where there are 2 J equal subdivisions in each direction. We solve the discrete obstacle 
problem on each level j with 8 = 1/3 so that the mesh parameter hj = (2 J — 1/3) -1 . 

We denote the energy norm on the j-th level by || • \\j. Let Uj be the numerical solution 
of the j-th level discrete obstacle problem and ej = UjU — Uj where IT, is the interpolation 
operator on the j— th level. We evaluate the error ||ej||j in the energy norm, and the error 
1 1 e j- 1 1 oo in the ioo norm, and compute the rates of convergence in these norms by 

P h = ha.{\\ej/2\\j/2/\\ej\\j)/hi{hj/2/hj) and (3^ = ln(||e i / 2 || 0O /||e i || 0O )/ln(/i i / 2 //ij)- 



The numerical results are presented in Table 4^ It is observed that the magnitude of the 
error in energy norm is O(h). 



3 


lejHj/IKHs 0h 


1 1 &j 1 1 oo /3oo 


1 


0.0000 xHT° 


0.0000 xio-° 


2 


1.2365 xlO" 1 


8.8312 xlO" 4 


3 


6.3226 xlO" 2 0.9094 


6.0088 xlO" 4 0.5221 


4 


2.5977 xl0- 2 1.2447 


8.8401 xlO" 5 2.6817 


5 


1.2159 xl0~ 2 1.0787 


2.4443 xl0~ 5 1.8267 


6 


5.9045 xl0~ 3 1.0343 


6.7946 xKT 6 1.8331 


7 


2.9125 xlO" 3 1.0157 


1.4775 xlO" 6 2.1929 


8 


1.4396 xlO" 3 1.0147 


8.8608 xKT 7 0.7363 



TABLE 4.1. Energy norm and norm errors for Example 1. 
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The exact coincidence set / for this example is the disc centered at (0,0) with radius Tq. 
Let Vj be the set of nodes on the j-th level corresponding to degrees of freedom involving 
pointwise evaluation of local basis functions in the interior of Q. Then we define the discrete 
coincidence set Ij by 

Ij = {P^ Vj ■ Uj(p) - ipip) < Halloo}- 



The discrete coincidence sets I 7 and I 8 are displayed in Figure |4.1[ where the radius of the 
circle in black is r . The convergence of the discrete coincidence sets is observed. 





Figure 4.1. Discrete coincidence set for Example 1 for level 7 (left) and level 8 (right). 
One of the advantages of the GFEM is that the local approximation space can be easily 



adjusted. In Table 4.2 we report the numerical results for the same problem but with Q 3 as 
the local approximation space. An 0(h 1,5 ) energy error is observed, which is due to the fact 
that the exact solution u is piecewise smooth. 



3 


ll e il j/l Mis 0h 


1 1 &j \ \ oo /3oo 


1 


1.4199 xl0~ 2 


1.0561 xlO" 4 


2 


6.1489 xlO" 2 -1.8589 


7.6281 xHT 4 -2.5078 


3 


1.8374 xlO- 2 1.6377 


9.4507 xnr 5 2.8313 


4 


5.8004 xHT 3 1.6134 


1.4396 xlO" 5 2.6330 


5 


2.3728 xl0~ 3 1.2702 


4.7114 xl0~ 6 1.5872 


6 


8.3768 xlO" 4 1.4908 


4.1685 xKT 7 3.4723 


7 


2.7675 xlO" 4 1.5918 


3.7129 xlO" 6 -3.1431 



TABLE 4.2. Energy norm and 
approximation space. 



norm errors for Example 1 with Q3 as local 



Remark 4.1. Note that the ioo norm errors fluctuate. This is likely due to the fact that the 
primal dual active set strategy is based on stopping conditions that are unrelated to the 
norm. 
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Example 2. In this example we take Q = (—0.5, 0.5) 2 , / = g = and ip(x) = 1 — 5|x| 2 + 
\x\ 4 . We solve the discrete obstacle problems using the same PU functions as in Example 1. 

Since the exact solution is not known, we take Cj = HjUj-i — Uj and compute the rates of 
convergence 0h an d Poo by 

h = ln(||e i /2||j/2/||ej|| i )/hi(/i j / 2 //ij) and (3^ = ln(||e i/2 ||oo/||ej W^)/ ln(h j/2 /hj). 
The results are presented in Table |4.3| 



3 


Il e jlli/ll%ll8 Ph 


.K'j || oc Poo 


1 


2.9288 xl0~° 


9.0040 x xlO" 1 


2 


5.9820 xlO" -0.9058 


5.3416 xx 10" 1 0.6622 


3 


1.2402 xl0~° 2.1333 


5.2357 xlO" 1 0.0271 


4 


6.5242 xlO" 1 0.8988 


2.5914 xlO" 2 4.2061 


5 


1.8496 xHT 1 1.7913 


1.7757 xl0~ 3 3.8091 


6 


8.9273 xl(T 2 1.0430 


4.4337 xHT 4 1.9867 


7 


4.4296 xl(T 2 1.0072 


1.1284 xl0~ 4 1.9667 


8 


2.2154 xKT 2 0.9977 


3.7776 xl0~ 5 1.5758 



TABLE 4.3. Energy norm and norm errors for Example 2. 



Since A 2 ip — f > in this example, the non-coincidence set is known to be connected [15] . 



This is confirmed by the discrete coincidence sets I7 and Is displayed in Figure 4^2 Note 
that the discrete coincidence sets have the correct symmetries: rotations by right angles and 
reflections across coordinates axes. 





■0£ 0.2 



Figure 4.2. Discrete coincidence set for Example 2 for level 7 (left) and level 8 (right). 
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j 


I e?l ?/l \ur\ 8 $h 

II J \ \J f II 1 1 <J ' 




1 


3.0796 xl0-° 


8.9960 xlO" 1 


2 


6.2833 xl(T -0.9044 


4.8507 xlO" 1 0.7834 


3 


1.0279 xlO" 2.4544 


4.3181 xlO" 1 0.1576 


4 


2.9125 xlO" 1 1.7646 


1.9025 xlO" 2 4.3689 


5 


1.4890 xlO- 1 0.9533 


1.6296 xlO" 3 3.4920 


6 


7.1583 xl0~ 2 1.0487 


4.8682 xlO" 4 1.7300 


7 


3.6108 xlO" 2 0.9836 


1.3055 xlO" 4 1.8916 


8 


1.8072 xlO" 2 0.9966 


3.1174 xlO" 5 2.0623 



TABLE 4.4. Energy norm and norm errors for Example 3. 



Example 3. In this example we take f2 = (—0.5, 0.5) 2 , / = g — and ip(x) = 1 — 5|x| 2 — 
\x\ A . We solve the discrete obstacle problems using the same PU functions as in Example 1. 
Numerical results are tabulated in Table 14.41 

The set-up for Example 3 is very similar to that of Example 2, except that now A 2 ip — f < 
and hence the interior of the coincidence set must be empty, otherwise the complementarity 
form of the variational inequality would be violated. This is confirmed by the discrete 



coincidence sets in Figure 4.3, which also possess the correct symmetries. 





0.2 0.2 



■0.4 -0.2 



0.2 0.4 



Figure 4.3. Discrete coincidence set for Example 3 for level 7 (left) and level 8 (right). 



Example 4. In this example we take Q to be the L-shaped domain (— 0.5, 0.5) 2 \[0, 0.5] 2 , 

We solve the discrete obstacle problems using 



(zi+0.25) 2 
0.2 2 



+ 



0.35 2 



/ = g = and ip(x) = 1 

a similar partition as described in Section 2A For this example, j is chosen so that it is the 
level where there are 2 J + 1 subdivisions in each direction, making hj = (2 J ' + 1 — 1/3) -1 . This 
allows us to insert an L-shaped element in the vicinity of the reentrant corner as described 
in Remark 12.31 
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From the numerical results in Table 4.5 we observe that (3h is approaching 0(h a ) where 
a = 0.544 is the index of elliptic regularity for the L-shaped domain, as predicted by Theorem 



3 


IK'Hi/IKHs 0h 


1 1 C j 1 1 oo fioo 


1 


4.4737 xl(T 


1.0000 xlO -0 


2 


6.9545 xl(T -0.7884 


5.9996 xlO" 1 0.9129 


3 


2.9079 xl(T 1.4086 


3.2598 xlO" 1 0.9854 


4 


1.8562 xlO -0 0.6864 


1.3853 xlO" 1 1.3086 


5 


6.9086 xlO" 1 1.4687 


4.0400 xlO" 2 1.8312 


6 


2.8930 xlO" 1 1.2747 


2.9381 xKT 2 0.4664 


7 


1.6919 xlO" 1 0.7797 


1.4457 xlO" 2 1.0308 


8 


1.0582 xlO" 1 0.6796 


6.9259 xlO" 3 1.0657 



TABLE 4.5. Energy norm and ioo norm errors for Example 4. 



Since A 2 ip — f = for this example, the non-coincidence set is connected [15], which is 



confirmed by Figure 4.4 



0.2 0.2 



-0.4 -0.2 0.2 



Figure 4.4. Discrete coincidence set for Example 4 for level 7 (left) and level 8 (right). 

Example 5. In this example we take Q to be the pentagon {x G (—0.5, 0.5) 2 : X\ + x<± < 
0.5}. We take / = g = and ip(x) = 1 — 9\x\ 2 . We solve the discrete obstacle problems using 
a similar partition as described in Section 2.1 For this example, j is chosen so that it is the 
level where there are 2 J + 1 subdivisions in each direction, making hj = (2 J + 1 — 1/3) _1 . This 



allows us to insert different types of elements near the obtuse vertices of fi, see Figure 4.5 



and Figure 4i3 The numerical results are reported in Table 4.6 

Since A 2 ip — f = in this example, the non-coincidence set is connected [15], which is 



confirmed by Figure 4.7 where the discrete coincidence sets also display the correct reflection 
symmetry. 
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Figure 4.5. Reference elements for the pentagonal domain. 




mm — v 



--?- T -•"-?- -T - 
T T • T T 

- -4 - f -- -f | 4- 

- -*- 4 ■ 
i 1 



• ~T" t " 
■ -1- I ■ 



• - 
I 1 

-T- -t - 
-1-4- 



T * 

. -4- 

* t 



/ rrrrrrfrrn v 



Figure 4.6. A partition of the pentagonal domain f2 for levels j = I (left) 
and j = 2 (right). The solid lines separate the different patches Qj,j = 1, ... ,9 
(left) and Qj,j = 1, ... ,25 (right). The dashed lines represent the extension of 
Qj by 5(h/2) on each side. This figure also shows the locations of the degrees 
of freedom. 
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j 


I e?l 7/I \ur\ 8 0h 

II J \ \J f II 1 1 a ' '» 


lle-ll 8 


1 


5.1600 xl(T 
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2 
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